  /* this header file yields r2, scaled and excluded for computing nonbonded interactions with various methods*/
  /* if excluded particles also needs processing (e.g. subtracting long ranged coulomb), define NO_SKIP_EXCL before include*/
  real rcut = p->rcut;
  real inv_rcut = 1. / rcut;
  real rsw = p->rsw;
  real scale14 = p->scale14;
  real coul_const = p->coul_const;
  vdw_param_t *vdw_param = p->vdw_param;
  real nonb_threshold = p->nonb_threshold;
  real rcut2 = rcut * rcut;
  real rsw2 = rsw * rsw;
  real inv_rcut2 = 1. / rcut2;
  real denom_switch = rcut2 - rsw2;
  real inv_denom_switch = 1 / (denom_switch * denom_switch * denom_switch);
  int order = p->order;
  double totalw = 0, totalu = 0;

  FOREACH_LOCAL_CELL(grid, ii, jj, kk, icell) {
    vec<real> oi, hi;
    cell2box(icell, &oi, &hi);
    real epsi[CELL_CAP], sigi[CELL_CAP], epsi14[CELL_CAP], sigi14[CELL_CAP];
    for (int i = 0; i < icell->natom; i++) {
      epsi[i] = vdw_param[icell->t[i]].eps;
      sigi[i] = vdw_param[icell->t[i]].sig;
      epsi14[i] = vdw_param[icell->t[i]].eps14;
      sigi14[i] = vdw_param[icell->t[i]].sig14;
    }
    FOREACH_NEIGHBOR_HALF_SHELL(grid, ii, jj, kk, dx, dy, dz, jcell) {
      int far = abs(dx) + abs(dy) + abs(dz) > 1;
      int jcell_id = jcell - grid->cells;
      int did = hdcell(grid->nn, dx, dy, dz);
      vec<real> dcell;
      vecsubv(dcell, jcell->basis, icell->basis);
      // real epsi[CELL_CAP], sigi[CELL_CAP], epsi14[CELL_CAP], sigi14[CELL_CAP];

      /*this is a marker that indicates the most recent i atom has excl/scal with index j*/
      int is_excl[CELL_CAP];
      int is_scal[CELL_CAP];

      for (int i = 0; i < icell->natom; i++) {
        is_excl[i] = -1;
        is_scal[i] = -1;
      }

      int ie = icell->first_excl_cell[did];
      int is = icell->first_scal_cell[did];
      for (int j = 0; j < jcell->natom; j++) {
        vec<real> xj;
        vecaddv(xj, dcell, jcell->x[j]);
        if (far && dist2_atom_box(&xj, &oi, &hi) > rcut2)
          continue;
        while (ie < icell->first_excl_cell[did + 1] && icell->excl_id[ie][1] <= j) {
          is_excl[icell->excl_id[ie][0]] = icell->excl_id[ie][1];
          ie++;
        }
        while (is < icell->first_scal_cell[did + 1] && icell->scal_id[is][1] <= j) {
          is_scal[icell->scal_id[is][0]] = icell->scal_id[is][1];
          is++;
        }
        real sigj = vdw_param[jcell->t[j]].sig;
        real epsj = vdw_param[jcell->t[j]].eps;
        real sigj14 = vdw_param[jcell->t[j]].sig14;
        real epsj14 = vdw_param[jcell->t[j]].eps14;
        int ni = icell->natom;
        if (dx == 0 && dy == 0 && dz == 0)
          ni = j;
        for (int i = 0; i < ni; i++) {
#ifndef NO_SKIP_EXCL
          if (is_excl[i] == j)
            continue;
#else
          int excluded = (is_excl[i] == j);
#endif
          int scaled = (is_scal[i] == j);
          vec<real> datom;
          vecsubv(datom, icell->x[i], xj);
          real r2 = vecnorm2(datom);
          if (r2 < rcut2) {